{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Building Chemical Reaction Networks (CRNs) Directly with BioCRNpyler\n",
    "\n",
    "**Overview:**  This tutorial shows how to use [BioCRNpyler](https://github.com/BuildACell/BioCRNPyler) to represent simple CRNs\n",
    "\n",
    "## What is a CRN?\n",
    "A CRN is a widely established model of chemistry and biochemistry.\n",
    "* A set of species $S$\n",
    "* A set of reactions $R$ interconvert species $I_r$ to $O_r$\n",
    "\n",
    "\\begin{align}\n",
    "\\\\\n",
    "I \\xrightarrow[]{\\rho(s)} O\n",
    "\\\\\n",
    "\\end{align}\n",
    "\n",
    " * $I$ and $O$ are multisets of species $S$. \n",
    " * $\\rho(s): S \\to \\mathbb{R}$ is a function that determines how fast the reaction occurs.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Import everything from biocrnpyler\n",
    "from biocrnpyler import *"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Combining Species and Reactions into a CRN\n",
    "\n",
    "The following code defines a species called 'S' made out of material 'material'. Species can also have attributes to help identify them. Note that Species with the same name, but different materials or attributes are considered different species in terms of the reactions they participate in.\n",
    "\n",
    "    S = Species('name', material_type = 'material', attributes = [])\n",
    "\n",
    "The collowing code produces a reaction R\n",
    "    \n",
    "    R = Reaction(Inputs, Outputs, k)\n",
    "\n",
    "here Inputs and Outputs must both be a list of Species. the parameter k is the rate constant of the reaction. By default, propensities in BioCRNpyler are massaction:\n",
    "\n",
    "### $\\rho(S) = k \\Pi_{s} s^{I_s}$\n",
    "\n",
    "Note: for stochastic models mass action propensities are $\\rho(S) = k \\Pi_{s} s!/(s - I_s)!$.\n",
    "\n",
    "Massaction reactions can be made reversible with the k_rev keyword:\n",
    "\n",
    "    R_reversible = Reaction(Inputs, Outputs, k, k_rev = krev)\n",
    "\n",
    "is the same as two reactions:\n",
    "\n",
    "    R = Reaction(Inputs, Outputs, k)\n",
    "    Rrev = Reaction(Outputs, Inputs, krev)\n",
    "\n",
    "\n",
    "Finally, a CRN can be made by combining species and reactions:\n",
    "\n",
    "    CRN = ChemicalReactionNetwork(species = species, reactions = reactions, initial_condition_dict = {})\n",
    "\n",
    "Here, initial_condition_dict is an optional dictionary to store the initial values of different species. \n",
    "\n",
    "    initial_condition_dict = {Species:value}\n",
    "\n",
    "Species without an initial condition will default to 0."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Species can be printed to show their string representation: m1_A_attribute m1_B m2_B D\n",
      "\n",
      "Reactions can be printed as well:\n",
      " m1[A(attribute)] --> 2m1[B] \n",
      " m1[B] --> m2[B]+D\n",
      "\n",
      "Directly printing a CRN shows the string representation of the species used in BioCRNpyler:\n",
      "Species = m1_A_attribute, m1_B, m2_B, D\n",
      "Reactions = [\n",
      "\tm1[A(attribute)] --> 2m1[B]\n",
      "\tm1[B] --> m2[B]+D\n",
      "]\n",
      "\n",
      "CRN.pretty_print(...) is a function that prints a more customizable version of the CRN, but doesn't show the proper string representation of species.\n",
      "Species (4) = {0. m1[A(attribute)] init_conc = 0, 1. m1[B] init_conc = 0, 2. m2[B] init_conc = 0, 3. D init_conc = 0}\n",
      "\n",
      "Reactions (2) = [\n",
      "0. m1[A(attribute)] --> 2m1[B]\n",
      " Kf=k_forward * m1_A_attribute\n",
      "  k_forward=3.0\n",
      "\n",
      "1. m1[B] --> m2[B]+D\n",
      " Kf=k_forward * m1_B\n",
      "  k_forward=1.4\n",
      "\n",
      "]\n"
     ]
    }
   ],
   "source": [
    "#Example: Model the CRN consisting of: A --> 2B, 2B <--> B + C where C has the same name as B but a new material\n",
    "A = Species(\"A\", material_type = \"m1\", attributes = [\"attribute\"])\n",
    "B = Species(\"B\", material_type = \"m1\")\n",
    "C = Species(\"B\", material_type = \"m2\")\n",
    "D = Species(\"D\")\n",
    "\n",
    "print(\"Species can be printed to show their string representation:\", A, B, C, D)\n",
    "\n",
    "#Reaction Rates\n",
    "k1 = 3.\n",
    "k2 = 1.4\n",
    "k2rev = 0.15\n",
    "\n",
    "#Reaciton Objects\n",
    "R1 = Reaction.from_massaction([A], [B, B], k_forward = k1)\n",
    "R2 = Reaction.from_massaction([B], [C, D], k_forward = k2)\n",
    "\n",
    "print(\"\\nReactions can be printed as well:\\n\", R1,\"\\n\", R2)\n",
    "\n",
    "#create an initial condition so A has a non-zero value\n",
    "initial_concentration_dict = {A:10}\n",
    "\n",
    "#Make a CRN\n",
    "CRN = ChemicalReactionNetwork(species = [A, B, C, D], reactions = [R1, R2], initial_concentration_dict = initial_concentration_dict)\n",
    "\n",
    "#CRNs can be printed in two different ways\n",
    "print(\"\\nDirectly printing a CRN shows the string representation of the species used in BioCRNpyler:\")\n",
    "print(CRN)\n",
    "\n",
    "print(\"\\nCRN.pretty_print(...) is a function that prints a more customizable version of the CRN, but doesn't show the proper string representation of species.\")\n",
    "print(CRN.pretty_print(show_materials = True, show_rates = True, show_attributes = True))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## CRNs can be saved as SBML and simulated\n",
    "\n",
    "To save a CRN as SBML:\n",
    "\n",
    "    CRN.write_sbml_file(\"file_name.xml\")\n",
    "\n",
    "To simulate a CRN with biosrape:\n",
    "\n",
    "    Results, Model = CRN_expression.simulate_with_bioscrape(timepoints, initial_condition_dict = x0)\n",
    "\n",
    "Where x0 is a dictionary: x0 = {species_name:initial_value}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Saving and simulating a CRN\n",
    "CRN.write_sbml_file(\"build_crns_directly.xml\")\n",
    "\n",
    "\n",
    "try:\n",
    "    import bioscrape\n",
    "    import numpy as np\n",
    "    import pylab as plt\n",
    "    import pandas as pd\n",
    "    \n",
    "    #Initial conditions can be set with a dictionary:\n",
    "    x0 = {A:120}\n",
    "    timepoints = np.linspace(0, 1, 100)#Timepoints to simulate over\n",
    "    \n",
    "    #This function can also take a filename keyword to save the file at the same time\n",
    "    R = CRN.simulate_with_bioscrape_via_sbml(timepoints = timepoints, initial_condition_dict = x0)\n",
    "\n",
    "    #Check to ensure simulation worked\n",
    "    #Results are in a Pandas Dictionary and can be accessed via string-names of species\n",
    "    plt.plot(R['time'], R[str(A)], label = \"A\")\n",
    "    plt.plot(R['time'], R[str(B)], label = \"B\")\n",
    "    plt.plot(R['time'], R[str(C)], \"--\", label = \"C\")\n",
    "    plt.plot(R['time'], R[str(D)],\":\", label = \"D\")\n",
    "    plt.legend()\n",
    "    \n",
    "except ModuleNotFoundError:\n",
    "    print(\"Plotting Modules not installed.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## ComplexSpecies and OrderedComplexSpecies\n",
    "\n",
    "When Species bind together to form complexes, it is recommended to use the Complex function:\n",
    "\n",
    "    Complex([list of species], ordered = True/False) \n",
    "\n",
    "This function returns the classes ComplexSpecies or OrderedComplexSpecies subclasses which contain information about the species inside of them. ComplexSpecies treats its internal species as an unordered multiset. OrderedComplexSpecies treats its internal species as an ordered list. It is recommended to always use the function Complex to create these types of Species for compatability reasons discussed in the OrderedPolymerSpecies example notebook.\n",
    "\n",
    "_Note: These objects do not automatically generate binding reactions. To do that, use the Component wrappers ChemicalComplex and OrderedChemicalComplex._"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "For ComplexSpecies, the order of the elements does not matter:\n",
      "C1=ComplexSpecies([A, B, B2, A])= complex_m1_A_2x_m1_B_m2_B_\n",
      "C2=ComplexSpecies([B, A, B2, A])= complex_m1_A_2x_m1_B_m2_B_\n",
      "C1==C2 ==> True\n",
      "\n",
      "For OrderedComplexSpecies, the Order Does Matter:\n",
      "C3=OrderedComplexSpecies([A, B, B2, A])= ordered_complex_m1_A_m1_B_m2_B_m1_A_\n",
      "C4=OrderedComplexSpecies([B, A, B2, A])= ordered_complex_m1_B_m1_A_m2_B_m1_A_\n",
      "C3==C4 ==> False\n",
      "\n",
      "ComplexSpecies (and the OrderedComplexSpecies and Multimers) are Species and can be used in reactions:\n",
      "Reaction.from_massaction([A, B, B2, A], [C1], k_forward = 10)=\n",
      "2m1[A]+m1[B]+m2[B] --> complex[2x_m1[A]:m1[B]:m2[B]]\n"
     ]
    }
   ],
   "source": [
    "A = Species(\"A\", material_type = \"m1\")\n",
    "B = Species(\"B\", material_type = \"m1\")\n",
    "B2 = Species(\"B\", material_type = \"m2\")\n",
    "\n",
    "print(\"\\nFor ComplexSpecies, the order of the elements does not matter:\")\n",
    "C1 = Complex([A, B, B2, A])\n",
    "C2 = Complex([B, A, B2, A])\n",
    "print(\"C1=ComplexSpecies([A, B, B2, A])=\", C1)\n",
    "print(\"C2=ComplexSpecies([B, A, B2, A])=\", C2)\n",
    "print(\"C1==C2 ==>\", C1==C2)\n",
    "\n",
    "print(\"\\nFor OrderedComplexSpecies, the Order Does Matter:\")\n",
    "C3 = Complex([A, B, B2, A], ordered = True)\n",
    "C4 = Complex([B, A, B2, A], ordered = True)\n",
    "print(\"C3=OrderedComplexSpecies([A, B, B2, A])=\", C3)\n",
    "print(\"C4=OrderedComplexSpecies([B, A, B2, A])=\", C4)\n",
    "print(\"C3==C4 ==>\", C3==C4)\n",
    "\n",
    "print(\"\\nComplexSpecies (and the OrderedComplexSpecies and Multimers) are Species and can be used in reactions:\")\n",
    "R = Reaction.from_massaction([A, B, B2, A], [C1], k_forward=10)\n",
    "print(\"Reaction.from_massaction([A, B, B2, A], [C1], k_forward = 10)=\")\n",
    "print(R)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Non-massaction propensities in BioCRNpyler\n",
    "By default, BioCRNpyler assumes that propensities are massaction with only one parameter, the rate constant $k_{forward}$ (and optional reverse rate $k_{reverse}$. However, additional propensity types are also supported. These reactions are created in two steps: first a Propensity of the appropriate type is made, then Reaction is made using that Propensity. These reactions are always created irreversibly. Examples are shown below:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## HillPositive: \n",
    "$\\rho(s) = k \\frac{s_1^n}{K^n+s_1^n}$\n",
    "\n",
    "Requried parameters: rate constant \"k\", offset \"K\", hill coefficient \"n\", hill species \"s1\"."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "m1[A] --> m1[B]\n",
      " Kf = k R^n / ( K^n + R^n )\n",
      "  k=1\n",
      "  K=5\n",
      "  n=2\n",
      "\n"
     ]
    }
   ],
   "source": [
    "#create the propensity\n",
    "R = Species(\"R\")\n",
    "hill_pos  = HillPositive(k=1, s1=R, K=5, n=2)\n",
    "\n",
    "#create the reaction\n",
    "r_hill_pos = Reaction([A], [B], propensity_type = hill_pos)\n",
    "\n",
    "#print the reaction\n",
    "print(r_hill_pos.pretty_print())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## HillNegative: \n",
    "$\\rho(s) = k \\frac{1}{K^n+s_1^n}$\n",
    "\n",
    "Requried parameters: rate constant \"k\", offset \"K\", hill coefficient \"n\", hill species \"s1\"."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "m1[A] --> m1[B]\n",
      " Kf = k R^n / ( K^n + R^n )\n",
      "  k=1\n",
      "  K=5\n",
      "  n=2\n",
      "\n"
     ]
    }
   ],
   "source": [
    "#create the propensity\n",
    "R = Species(\"R\")\n",
    "hill_neg  = HillPositive(k=1, s1=R, K=5, n=2)\n",
    "\n",
    "#create the reaction\n",
    "r_hill_neg = Reaction([A], [B], propensity_type = hill_neg)\n",
    "\n",
    "#print the reaction\n",
    "print(r_hill_neg.pretty_print())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## ProportionalHillPositive: \n",
    "$\\rho(s, d) = k d \\frac{s_1^n}{K^n + s_1^n}$\n",
    "\n",
    "Requried parameters: rate constant \"k\", offset \"K\", hill coefficient \"n\", hill species \"s1\", proportional species \"d\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "m1[A] --> m1[B]\n",
      " Kf = k D R^n / ( K^n + R^n )\n",
      "  k=1\n",
      "  K=5\n",
      "  n=2\n",
      "\n"
     ]
    }
   ],
   "source": [
    "#create the propensity\n",
    "R = Species(\"R\")\n",
    "D = Species(\"D\")\n",
    "prop_hill_pos  = ProportionalHillPositive(k=1, s1=R, K=5, n=2, d = D)\n",
    "\n",
    "#create the reaction\n",
    "r_prop_hill_pos = Reaction([A], [B], propensity_type = prop_hill_pos)\n",
    "\n",
    "#print the reaction\n",
    "print(r_prop_hill_pos.pretty_print())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## ProportionalHillNegative: \n",
    "$\\rho(s, d) = k d \\frac{1}{K^n + s_1^n}$\n",
    "\n",
    "Requried parameters: rate constant \"k\", offset \"K\", hill coefficient \"n\", hill species \"s1\", proportional species \"d\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "m1[A] --> m1[B]\n",
      " Kf = k D / ( K^n + R^2 )\n",
      "  k=1\n",
      "  K=5\n",
      "  n=2\n",
      "\n"
     ]
    }
   ],
   "source": [
    "#create the propensity\n",
    "R = Species(\"R\")\n",
    "D = Species(\"D\")\n",
    "prop_hill_neg  = ProportionalHillNegative(k=1, s1=R, K=5, n=2, d = D)\n",
    "\n",
    "#create the reaction\n",
    "r_prop_hill_neg = Reaction([A], [B], propensity_type = prop_hill_neg)\n",
    "\n",
    "#print the reaction\n",
    "print(r_prop_hill_neg.pretty_print())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## General Propensity: \n",
    "$\\rho(s) = $ function of your choice\n",
    "\n",
    "For general propensities, the function must be written out as a string with all species and parameters declared."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "m1[A]+m1[B] --> complex[2x_m1[A]:m1[B]:m2[B]]\n",
      "k1*2 - k2/S^2\n",
      "  k1=1.11\n",
      "  k2=2.22\n",
      "\n"
     ]
    }
   ],
   "source": [
    "#create species\n",
    "# create some parameters - note that parameters will be discussed in the next lecture\n",
    "k1 = ParameterEntry(\"k1\", 1.11)\n",
    "k2 = ParameterEntry(\"k2\", 2.22)\n",
    "S = Species(\"S\")\n",
    "\n",
    "#type the string as a rate then declare teh species and parameters\n",
    "general = GeneralPropensity(f'k1*2 - k2/{S}^2', propensity_species=[S], propensity_parameters=[k1, k2])\n",
    "\n",
    "r_general = Reaction([A, B], [C1], propensity_type = general)\n",
    "\n",
    "print(r_general.pretty_print())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}